#include <cstring>
#include <iostream>
#include <cmath>
#include <boost/math/distributions/students_t.hpp>
extern "C" {
    //#include <cblas.h>
    //#include <clapack.h>
#include <mkl/mkl.h>
} using namespace std;
/*
 * nx is the number of instances
 * mx is the number conditions
 * X[mx,nx]
 * Y[my,nx]
 * param[mx,my]
 * my is the number targets
 */

int regression(float* X, int mx, int nx, float* y, int my, float* param)
{
    boost::math::students_t dist(nx - 2);    //nx-2 is the degree of freedom
    float* XT = new float[mx * nx];
    float* XTX = new float[mx * mx];
    float* PREFIX = new float[mx * nx];
    //param = new float[my * mx]; This should be allocated before calling
    memcpy(XT, X, sizeof(float) * mx * nx);
    //C := alpha*op(A)*op(B) + beta*C
    float alpha = 1.0;
    float beta = 0;
    //XTX = X'(size mx*nx)*X(size nx*mx)
    //void cblas_sgemm (const CBLAS_LAYOUT Layout, const CBLAS_TRANSPOSE transa, const CBLAS_TRANSPOSE transb, const MKL_INT m, const MKL_INT n, const MKL_INT k, const float alpha, const float *a, const MKL_INT lda, const float *b, const MKL_INT ldb, const float beta, float *c, const MKL_INT ldc);
    cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
        mx, mx, nx, alpha, XT, mx, X, mx, beta, XTX, mx);
    //(XX')^-1  size mx * mx
    // - Compute the LU factorization of a M by N matrix A
    int* ipiv = new int[mx];
    LAPACKE_sgetrf(LAPACK_ROW_MAJOR, mx, mx, XTX, mx, ipiv);
    //When using atlas and openblas
    //clapack_sgetrf(CblasRowMajor, mx, mx, XTX, mx, ipiv);
    // - Generate inverse of the matrix given its LU decompsotion
    LAPACKE_sgetri(LAPACK_ROW_MAJOR, mx, XTX, mx, ipiv);
    //When using atlas and openblas
    //clapack_sgetri(CblasRowMajor, mx, XTX, mx, ipiv);
    //Now XTX contains (XX')^-1 size mx * mx
    // calculate (XX')^-1*X' size mx * nx
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, mx, nx, mx,
        alpha, XTX, mx, XT, mx, beta, PREFIX, nx);
    //Now PREFIX contains (XX')^-1*X' size mx * nx
    //begin calculate the parameters;
    float* ry = new float[nx];
    float* tmp = param;
    for (int i = 0; i < my; i++)
    {
        param += (i * (mx + 1 + mx - 1));
        y += (i * nx);
        //param=PREFIX*Y
        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, mx, 1, nx,
            alpha, PREFIX, nx, y, 1, beta, param, 1);
        //Calculate variance
        //First, calculate the predicted y value
        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nx, 1, mx,
            alpha, X, mx, param, 1, beta, ry, 1);
        float se = 0;
        for (int j = 0; j < nx; j++)
        {
            se += (ry[j] - y[j]) * (ry[j] - y[j]);
        }
        param[mx] = sqrt(se / nx);    //the standard sigma of error
    }
    param = tmp;
    for (int k = 0; k < my; k++)
    {
        param += (k * (mx + 1 + mx - 1));
        for (int j = 0; j < mx - 1; ++j)
        {
            float mean = 0;
            for (int i = 0; i < nx; ++i)
            {
                mean += X[i * mx + j + 1];
            }
            mean = mean / nx;
            float tvalue = 0;
            for (int i = 0; i < nx; ++i)
            {
                tvalue +=
                    (X[i * mx + j + 1] - mean) * (X[i * mx + j + 1] -
                        mean);
            }
            tvalue = fabs(param[mx] / sqrt(tvalue));
            tvalue = param[j + 1] / tvalue;
            tvalue = 1.0 - boost::math::cdf(dist, tvalue);    //+boost::math::cdf(dist, -tvalue);
            if (tvalue>=0.5) tvalue=1.0-tvalue;
            param[mx + j + 1] = 2 * tvalue;
        }
    }
    //now param contains the necessary parameters;
    delete[]XT;
    delete[]ipiv;
    delete[]XTX;
    delete[]PREFIX;
    delete[]ry;
    return 0;
}
